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Introduction 


Many methods of applied mathematics are based on the use of polynomial 
approximations, and in this unit we shall consider how some of these methods are 
derived and used to solve practical problems. In the next unit we shall make use of 
many of these results when solving differential equations numerically. 


We corisider two ways in which we might use a polynomial to approximate a 

function f(x). If we know the value of the function and its first n derivatives at 

some point x = a, say, as in Figure 1, then we can determine the Taylor In this unit the word point 

polynomial of degree n about « which approximates f(x) near the point x = a. means a real number in the 
domain of the function. 

In many problems the values of the derivatives are not known but the values of f 

are given at n + 1 points xo, X;, X2, ... , X, a8 in Figure 2. In this case we can 

determine that polynomial of degree n which takes these values of f at 


Xo,X1,-++X,- Such a polynomial is called an interpolating polynomial. 


J y 
Sia) 
fx) 
| Sx) x oJ 
mae or Sx) Ts = ee 
: ex f(x) 
$ x fx) 
s(x) x 
b a 7. T r T t = 
Figure 1. For Taylor polynomials xe Bai ks ae 
the data is given at x = , Figure 2. For interpolating polynomials the data is given at n + 1 points. 


Taylor polynomials form the basis for two important methods in this unit: a 
method for finding roots of equations and a method for finding the approximate 
solutions of differential equations. Further, in the television programme we look at 
how Taylor polynomials can be used to derive approximate’ relationships between 
variables where it is impossible to find an exact relationship. 


In the final section of the unit we use interpolating polynomials to derive methods 
for determining the approximate values of integrals, such as the trapezoidal 
method and Simpson's method. 


Study guide 


Section 1 of this unit revises the determination of Taylor polynomials and is 
central to the rest of the unit. Subsection 2.1 should be studied before Section 3, 
while Section 4 can be studied at any time after Section 1. Subsection 2.2, the tape 
subsection, can be studied at any time after Section 1, and is an important 
introduction to the next unit. Section 5 contains some extra exercises which you 
may attempt if you need more practice on the techniques introduced in this unit. 


Before watching the television programme you should have studied Subsection 2.1, 
and Section 3 up to Subsection 3.3. The television programme concerns the study 
of freely hanging ropes and cables. 


1 Basic theory 


1.1 Polynomials 
A polynomial is a function, p(x) say, which can be expressed in the form 
P(x) = ag + ax + apx? +--+ +.4,x" 


where n is a non-negative integer and do, a,, ... , @, are real numbers called the 
coefficients of the polynomial. 
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If nis the highest power of x for which the coefficient a, is non-zero, p(x) is said to 
be a polynomial of degree n. Hence the polynomial 

1+ 7x + 4x? 
is a polynomial of degree 3, 


Polynomials of low degree have acquired special names. 


These definitions are also 
given in Subsection 3.1 of the 
Handbook. 


Thus the polynomial 1 +7x + 4x° is referred to as a cubic polynomial or simply a 
cubic, 


The advantage of using polynomials is that they are very easy to manipulate 
algebraically. For example, if p(x) is given by 


p(x) = 1 +3x + 4x? + x3 

and q(x) is another polynomial given by 
q(x) = 3+ 2x —x?, 

then we can, for example, add the two polynomials to get another polynomial 
p(x) + q(x) = 4 + Sx + 3x? +x. 

Polynomials are also very easy to differentiate or integrate. To differentiate 
P(x) = 1 + 3x + 4x? + x°, 

we simply differentiate each term separately to obtain 


a(x) = 3+ 8x + 3x7, 
while integrating p(x) gives 
foc ae = fo + 3x + 4x? + x*)dx 


= 32 43 Lr 
=C+x+5x +e te 


where C is an arbitrary constant. 


Example 1 
The function tan x may be approximated near x = 0 by the polynomial 
1 2 p(x) is the Sth Taylor 
P(x) =x + 3x? + polynomial about 0 for tan x, 
3 "is" Taylor polynomials will be 


Use this approximation to estimate discussed in Subsection 1.2. 


5 
f tan x dx. 
0 


0.5 
j tanxdx = li§ frye+e x) ax 
-[} +7 nw +] 


= 0.130556 to six decimal places. 


Solution 
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(The true solution, using the integration tables in Subsection 7.1 of the Handbook, 


eve 0.5 5 
f tanxax =| —tog.cosx | 
0 


= 0.130584 to six decimal places. 
Hence our approximation method gives a result which is correct to four decimal 
places.) 


When evaluating a polynomial for a particular value of x, the most efficient 
method to use is called nested multiplication. This method will also usually reduce 
rounding errors in the calculations. To see how the method works we consider the 
polynomial 

P(x) = 1 + 3x + 7x? + 4x3. 
If, for convenience, we reverse the order of the terms in the polynomial so that we 
write 

P(x) = 4x3 + 7x? + 3x + 1, 


then the nested multiplication method can be used to evaluate p(x) as 
P(x) = ((4x + 7)x + 3)x + 1 


The steps in the calculation are as follows. 


i iter 4 into tl tor. 
0) Enter 4 into the calculator. Titi aac 


(ii) Compute 4x + 7. sequence of operations, given 
iii i that x is stored in the memory 
kp sge te ESS and can be obtained using the 
(iv) Compute ((4x + 7)x +3) x x +1. RCL key, is 
This sequence of operations can be expressed using a recurrence relation as HERA 
Mo = 4 F] [Req] FB 
Uy = Xo + 7 x] iRCLI 1 Ii 
U2 = xu, +3 
us = xu) + 1, 
and us gives the value of the polynomial at x. 
Example 2 
Evaluate p(x) = 4x° + 7x? + 3x + 1 at x = 1.27. 
Solution 
Using the recurrence relation, we have 
Up =4 


uy = 1.27uy + 7 = 12.08 

uz = 1.27u, + 3 = 18.3416 

U3 = 1,.27u, + 1 = 24.293 832. 
Notice that this method of evaluating the cubic requires only 3 multiplications and 
3 additions. 


The following procedure box gives the recurrence relation for evaluating a 
polynomial of degree n. 


Nested multiplication 
To evaluate 
P(X) = yx" + Gy_yx"~! +--+ + ax + ag 
at a particular value of x, the nested multiplication method uses 


the recurrence relation 
Up 41 = XU, + Oy—p—4 


where Up = a, 
and u, gives the value of p(x). 
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Exercise 1 
This exercise concerns the polynomial 
P(x) = 2x* + 3x9 + x? - 7x42. 
(i) What is the derivative of p? 
(ii) What is the integral of p? 
(iii) Evaluate p(x) for x = 0.29 using nested multiplication. 
fiv) Evaluate the derivative of p at x = 0.77. 
[Solution on p. 50] 


1.2. Taylor polynomials 

There are many functions such as sin x, e*, log, x and so on which we have to 
evaluate for a given value of x when carrying out computer or calculator 
computations. The traditional method, which you may have used, was to look up 
the required value of sin x in a table of sines. It is however not feasible for 
computers or calculators to store tables of numbers, since it would use up a 
relatively large amount of computer time to locate the required value. What 
happens in practice is that we approximate functions like the sine function by 
polynomials, so that for a given value of x we approximate sin x by the value of 
the polynomial. 


There are several types of polynomial approximation, but traditionally the most 
popular are the Taylor polynomials. Suppose that we want to determine a 
polynomial p(x) of degree n which approximates some given function f(x) near 
x = 0. We write down p(x) as 
P(X) = dp + ayx + apx? + +++ + a,x", 

To derive the Taylor polynomial, we require p to be such that 

PO) =f(0), 

pO) =f'0) 

p"(0) =f"(0), 


p(0) =f"). 
This uniquely determines p, since 
P(O) =a, 80 that a = f(0); 
p'(0)=a,, so that a, = f'(0); 


‘x i) 
p%0)= ka, so that a, =/ 1 
Hence 
2 er 
P(x) = $00) + x f'0) + FIO) + += $0) a) 


This polynomial will normally give very good approximations to f near x = 0. 
(However, there are some functions, such as |x|, which cannot be approximated in 
this way since their derivatives do not exist at x = 0.) The above polynomial is 
known as the nth Taylor polynomial about 0 for the function f- 


Here are two examples of Taylor polynomials, 


Example 3 
Determine the 3rd Taylor polynomial about 0 for the function 


Sxyae 


This derivation is also given in 
the Taylor Series units of 
MIOI and MS283. 


You have already used one 
Taylor polynomial in Section 
3 of Unit 7. See also Section 4 
of the Preparatory Booklet, in 
which the alternative name 
nth Taylor approximation is 
used. 
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Solution 
Differentiating this function is easy, since 
d 

Ror 

so that 
Se)=e,  f"()=F and f"(x)=e% 

At x = 0, e* = 1, so that 
£0) =f'O) =f") =F"(0) = 1. 

Hence the 3rd Taylor polynomial about 0 approximating f(x) = e* is obtained by 

substituting into Equation (1), giving 

x 3 

2 3r 

In Figure 1 we have plotted the graphs of e* and p(x). Note that the two graphs 

are indistinguishable near x = 0. 


p(x) =14 x4 


yoltxe dregs? 


Figure 1. Graphs of and p(x) = 1 +x + 3x2 + 22% 


Example 4 
Consider the function 
F(x) = (1+ xV, 

where n is a positive integer. This is a polynomial of degree n and we would 
normally use the Binomial Theorem to obtain the expansion. However, we shall The Binomial psc 

ii i i it i iven in M101 A 
obtain the same polynomial expansion by looking at Taylor polynomials. a "i tn ectises tot ie 
The derivatives of f(x) are Preparatory Booklet. 


Sx) =n(l + xP" 
L(x) = n(n — 1) + xP? 


fx) =n! 
All higher derivatives are zero. 
At x =0, (1 + x)* = 1 for any value of k, so that 
SO=1, fO=n, f"O)=n(n—-1), ..., FO)=n! 
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Hence, using Equation (1), we have the nth Taylor polynomial about 0 for 
f(x) = (1 + x)" as 


a Ya 4 n(n — ae 2). fe 


Here p(x) is exact because (1 + x)" is a polynomial of degree n, and we can write 


P(x) =1+nx+ eo 


Qs at 4 me MODs , MO DH 2) 5 
2 3! 

In an exercise at the end of this subsection you will be asked to determine a 

Taylor polynomial approximation for (1 + x)"/?. 


at 


Very frequently we need approximations near some point x = where a is not 
necessarily equal to zero. Suppose that we look for a polynomial of the form 


P(X) = do + ay(x — a) + ay(x — a)? + --> +.a,(x — a)" 
such that 
P(x) =f (a) 
P'(x) =f'(a) 


ya) =f(a). 


Such a polynomial is called the nth Taylor polynomial about « for the function /. 
Again it is not difficult to deduce that the coefficients of this polynomial are 


“thy, 
a= a 


Thus we may define these Taylor polynomials as follows. 


k=0,1,...,m 


The nth Taylor polynomial about « approximating the function 
J (x) is defined by 


(x - 
2! 


2 = 
Plex) = f(a) + («= a) f(a) + Spray $4 FM poy 


To evaluate polynomials of the form 


P(x) = do + ay(x — x) + a(x — a)? + ++ +.a,(x — a), 
the nested multiplication procedure can be adapted as 
Une = (X= a)lly + dy—p—a 
where up = a, and u, gives the value of p(x). Thus polynomials in this form are as When evaluating them, do not 


easy to evaluate as are ordinary polynomials. pops h ease ah as 
Example 5 
Determine the 3rd Taylor polynomial about 1 approximating the function 

fQx)=e. 


Evaluate this polynomial at x = 1.55. 
Solution 
For n= 0, 1, 2... we have f(x) = e* and f™(1) =e. Thus, 
P(x) =e + (x — le + Hx — 1'e + Hx — 1Pe 
=e{1 + (x—1) + Hx - 17 + Mx — 13}. 


This approximation is compared with e* in Figure 2 (overleaf). In this case the 
graph of the polynomial is indistinguishable from that of e* near x = 1. 
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yy yre 


yre{l+(x- I+ de IP + Ele 1} 


Figure 2. Graphs of p(x) = e{1 + (x — 1) + 4(x — 1)? + Hx — 1)*} and e*. 


To evaluate p(x) at x = 1.55, compute (x — 1) = 0.55, Using nested 
multiplication to evaluate 1 + (x — 1) + (x — 1)? + k(x + 1), we have 
Uo = 0.166 666 67 


Uy = 0.559 + 0.5 = 0.591 666 67 

Uz = 0.55u, + 1 = 1,3254167 

uz = 0.55u2 + 1 = 1.7289792. 
Hence p(1.55) = eus = 4.699 85 to five decimal places. This is not a bad 
approximation to the correct answer, which is 4.711 47 to five decimal places. 


Exercise 2 


Determine the 2nd Taylor polynomial about 0 approximating f(x) = (1 + x)'/?. Evaluate 
this polynomial at x = 0.44 and compare your answer with the true solution. 


Exercise 3 
Show that the nth Taylor polynomial about 0 approximating f (x) = log. (1 + x) is 


mer 
epi cee ad 


px) =x— $2 + bt — - 


[Solutions to Exercises 2 and 3 on p. 50). 


1.3 Error bounds for Taylor polynomials 


Wherever we use these approximations it is wise to obtain some estimate of the 
error in the approximation at a given value of x. That is, we would like to know 
something about the magnitude of the error function ¢(x) defined by 


a(x) = p(x) — f(x). (2) 


For example, the error function for the third Taylor polynomial about 0 
approximating e*, given in Example 3, is 


a(x) =1l+x+h74+¢b0-e. 
(x) has been plotted in Figure 3. 


On my calculator I stored 


0.55 in the memory using the 
RCL key to recall it. Here is 


my key sequence 


(Gz) &) 


RCL! 


+E] 


we 


=] Bl 


RCL 


+] [1 


HE 
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y= 0.187x4 


\ 

\ 
\ 
\ 

\ 
1 
4] 


y= -0.187x 


Figure 3. Graph of e(x) = 1 +x + x2 + 8 — et 


For comparison, the graphs of the functions +0.187 x* are also shown. (The 
reason for choosing these particular functions will become clear when you read 
Example 6 on pp. 12-13.) 


The behaviour of «(x) near x = 0 in the above example should be expected because 
we have required that the first three derivatives of p(x) and e* should agree at 

x = 0. Hence (0) = e'(0) = e"(0) = e’"(0) = 0, and so the first three Taylor 
polynomial approximations to e(x) are all zero. 


Suppose now that we have obtained the nth Taylor polynomial P(x) about « 
approximating some given function f(x), and we wish to use this polynomial p to 
approximate fat a particular value of x. We need some idea of what the 
magnitude of the error is at this value of x. The theorem that follows gives a 
method of doing this. It is usually stated not in terms of the error function e(x) 
but in terms of a function called the remainder, defined as the function r(x) in the 
formula 


L(x) = plx) + r(x). 


Comparing with (2), we see that the remainder and the error function are related 
by 


r(x) = —e(x). (3) 


Taylor’s Theorem 

Let f and its first n + 1 derivatives exist and be continuous at all The requirement that the 
points between and including x and x. Then f(x) can be expressed n+ 1 derivatives of f(x) _ 

as should exist excludes functions 
such as |x| or x"!? over an 
(Co interval which includes 0. 


2 
5 Playas 


L(x) =f (a) + (x — a) f(a) + 


oo FED pm) + 100) 


where r(x), the remainder, satisfies 


r(x) = = oF is Yc.) T have used c, here rather 


(n+ 1)! than c to indicate that c, 


dey the value of x. 
for some number c, between « and x. pends on the value of x. 


The expression (4) is often called the nth Taylor expansion of f (x) about « with 
remainder. The infinite sum 


— a)? = 
Fes) =F a) + — af) + SS pray ¢ + FM pay --- 5) 


2 MST204 18.1 


is called the Taylor expansion (or Taylor series) of f(x) about «. Its significance is 
that, if we know for some particular value of x that the expression for the 
remainder, 
x (-oyt! “(n+ Ly 

r(x) = “eri (cx), 
approaches zero as n becomes large, then we can conclude that (for that value of 
x) the sum of the terms in the Taylor series converges to f(x). Although we will 
not be going into the question of convergence, the use of the Taylor series 
expression for a function is so widespread in applied mathematics that it is 
important to remember expression (5). 


Rather than prove Taylor's Theorem I will just show what happens in the 
particular case when n = 0. In this case Taylor’s Theorem states that if f and f’ 
exist and are continuous for all points between and including « and x, then 


L(x) =f (a) + (x — a) f (cx) (6) 
where c, lies between « and x. 
Rearranging Equation (6) gives 
f(x) — f(a) 


Ea =e 
Graphically this can be interpreted as saying that between « and x there must be a 
point c, for which the slope of the graph (f’(c,)) is the same as the slope of the 
(dashed) line joining the points (a, f(x)) and (x, f(x)) in Figure 4, 


line is 


pal 


Figure 4. Graph illustrating Taylor's Theorem for n = 0. 


This special case with n = 0 has an alternative name. It is often called the Mean 
Value Theorem for derivatives. 


Since the remainder and the error are equal and opposite, we can use Taylor's 
Theorem to estimate the error function e(x). The following example shows how 
this can be done. 


Example 6 

(0) Verify that for all c satisfying 0 < c < 1.5 we have |e| < e'* = 4.4817... . 

(ii) Deduce that for all x satisfying 0 < x < 1.5 the error of the 3rd Taylor 
polynomial about 0 for e satisfies 


le(x)| << Ax*, — where A = ae 


Solution 
(i) Since e* is always positive we have |e‘| = e*, and since c < 1.5 and the 
function is increasing we have e < e!-> = 4.4817... . 


It follows that |e| < 4.4817... as required. 
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(ii) Taylor’s Theorem, with f(x) = e*, gives 
@=14+x + $x? + be + r(x), where 
r(x) = dex* fX(c,) 
= tex* exp(c,). (7) Remember that exp (...) 


means e to th ices 
Since 0 < x <5, it follows that 0 <c, < 1.5 and hence, by part (i), that pee cr om 
lexp(c.)| <e"*. 


Using (7) we deduce that 


Ir(x)| < zex* eS 
= Ax‘, where A = ye! = 0.187. 


Since (x) = —r(x) (Equation (3)), we have |e(x)| = |r(x)|, and hence 
le(x)| < Ax* 
as required, 
This result can also be written 
—Ax* < e(x) < Ax*. 


Figure 3 (page 11) shows that these inequalities are indeed satisfied for 
0 <x < LS, as the theorem states. 


The procedure for determining the error bound used in Example 6 is formalized in 
the following box. 


Error bounds for Taylor polynomials 
1. Let p(x) be the nth Taylor polynomial about « approximating 
a given function f(x) which satisfies the conditions of Taylor's 
Theorem. We want to find a bound on the error function 
#(x) = p(x) — f(x) 


for all x satisfying a < x <b, where a, b are given numbers 
such that a <a and a <b. 


2. Determine a value for M such that 
Ife" %(c)| <M 
for all c satisfying a<c <b. 


3. Taylor's Theorem gives 


= ae 
- mr FO MCs) 


where c, lies between and x, and it follows since a < c, <b, 
that 


e(x) = 


Ix — a"?! 
(n+ 1)! 


This is the required error bound. 


le(x)| < M. 


Exercise 4 


(i) Verify that for 0 <c < 0.5 we have | +c] < 1.5. 


(ii) Compute an error bound for the 3rd Taylor polynomial about 0 approximating 
(1 + x)$, valid for all x satisfying 0 < x < 0.5. (Note that you are not asked to 
compute the 3rd Taylor polynomial itself.) 


Exercise 5 
(i) Verify that |sin cl <1 for all c. 


(ii) Determine the Sth Taylor polynomial about 0 approximating sin x. Compute an 
error bound for this approximation at x = I, using Taylor's Theorem. 


[Solutions to Exercises 4 and 5 on p. 50] 
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Summary of Section 1 
1. A polynomial is a function, p(x) say, which can be expressed in the form 
P(x) = ag + aX + a2Xx? + +> + ayx" 
where n is a non-negative integer and do, a;, ... , 4, are real numbers called the 
coefficients of the polynomial. If a, is non-zero, p is called a polynomial of degree n. 
2. To evaluate p(x), written as 
P(X) = yx" + dy 1X") + ++> + 1X + do, 


for a particular value’of x, the nested multiplication method uses the recurrence 
relation 


Up y= XM, + Oy —p—t 


where to = a, 
and u, gives the value of p(x). 


If p(x) is of the form 

P(x) = ay(x — a)" + Gq a(x — af! + = + a(x — 2) + ao 
then the recurrence relation is 

Uy = (X — Olly + Appa 
where, as before, uo = a, and u, gives the value of p(x). 


3. The nth Taylor polynomial about « approximating / is given by 


—ay =- 
phx) = F(a) + («= af (0) + SP pra) + ¢ R= AP pay 


4. Taylor’s Theorem. Let f(x) and its first n + 1 derivatives exist and be 
continuous at all points between and including « and x. Then f(x) can be 
expressed as 


= 2 
$e) = fle) + & — af (a) + SPX pr(a) + = 


= aP peng 4g R= pensry 
+ f(a) + (+! LON), 
—ayt! 
where c, lies between « and x and the expression at + 1(c,) is the 
remainder. 3 
Hence the error function ¢(x) can be written as 
a(x) = p(x) — f(x) 
= aa int ty 
res SO" MCs). 
If a number M can be found such that 
User Mc <M 
for all c satisfying a < c < b, then it follows that 
dol il 
ie Lererar aa 


This is an upper bound on the error. 
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2 Two applications of Taylor polynomials 


2.1 Roots of equations: the Newton-Raphson method 


In the course of solving mathematical problems it is often necessary to find the 
values of x for which the equation 


f(x) =0 
is satisfied, where f is some known function. I shall refer to them as the roots of They can also be called the 
the equation f(x) = 0. solutions of the equation. You 
: have already met these words, 
One way of determining the roots of an equation is to sketch the graph of used interchangeably, in 
oe Unit 5, Complex numbers. 
y=IX 


and to determine those values x for which the curve cuts the x-axis. 
For example, to find the roots of the equation 

f(x) = 2x? — 5x? — 5x +10=0 
we could sketch the graph as in Figure 1. 


Curve sketching is reviewed in 
Subsection 5,3 of the 


0 


at Handbook. 
=) 
y(-2)=-16 if 

“lo 
Figure 1, Sketch of y = 2x? — Sx? — 5x + 10. 
Looking at Figure 1 we can see that there are three values of x at which the graph 
cuts the x-axis. These are approximately —1.5, 1.2 and 2.8. 
In very rare cases the graph may only touch the x-axis but not cross it, so that the ‘ 
graph may be as in Figure 2. Figure 2. Graph of a function 

which touches but does not 
It can be difficult to determine roots in this instance, since with the naked eye we cross the x-axis. 
cannot detect whether there is one root, two roots or no roots as in the following 
close-ups. 
{a) one root (b) two roots (c) no roots 


Figure 3. Possible close-up views of the graph in Figure 2. 
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However in this unit we will assume that the graph is of the type given in Figure 1. 


Another method of finding roots, which you have met before, is the formula See Subsection 3.1 of the 
method for quadratic equations such as Handbook. 
x? —13x+7=0. 


Here are some equations for which the roots are not so easy to compute. 
(i) x—2sinx=0 

(ii) x? + 13.5x? + 40x + 16.67 =0 

(ii) xe*-0.1=0 


What we want is a method which could be used to solve any of these problems. 
The method we will use is an application of Taylor polynomials and is called the 


Newton-Raphson method. 

The Newton-Raphson method 

This method determines the roots of an equation using a recurrence relation. We ‘The Newton-Raphson method 

begin by supposing that there is a root of the equation f(x) = 0 near some point was introduced in the Taylor 

Xo. Now the first Taylor polynomial about xo approximating f(x) is given by eae of M101 and 
P(x) =f (Xo) + (x — Xo) f'(X0). () 


This approximation is often called the tangent approximation to f(x) at Xo because 
it is the equation of a straight line which is tangential to the curve y = f (x) at Xo, 
as in Figure 4, 


4 


ie) the tangent 
approximation 
al x, 


Px) = S09) + (= Xo) SCX9) 


av 


Figure 4. The tangent approximation to f(x) at xo. 


In Figure 4, Xo is not a particularly good approximation to the root of f(x) = 0, 
but we can use the tangent approximation at x» to improve this. The graph of the 
tangent approximation cuts the x-axis at the point x, such that 

P(x1) = 0. 
Since p is an approximation to f, the root, x,, of the equation p(x) = 0 is likely to 
be a close approximation to the root of f(x) = 0. 
To calculate x; we use Equation (1), which tells us that p(x,;) = 0 when 

f(%0) + (%1 — Xo)f (xo) = 0. 


Hence rearranging gives 
‘f(%o) A slightly different derivation 
ar (2) of Equation (2) will be given 
f'(%o) in the television programme. 
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Example 1 
From Figure 1 we know that there is a root of the equation 
F(x) = 2x? — 5x? — 5x +10=0 Remember to use nested 
multiplication to evaluate f(x), 
near Xo = —1.5. 


Now f (—1.5) = —0.5 
and f’(x) = 6x* — 10x — 5, giving 


f'(-1.5) = 23.5. 
Using Equation (2) with x9 = —1.5 gives 
i OES) ogo OS. 
x,=-15 f'(15) > 1S 4555 = —1.478723 4. 


We can see that x, is probably a much better approximation to the root than Xo, 
since f(—1.478 7234) = —0.006 318 39 and this is considerably smaller in 
magnitude than f(—1.5) = —0.5. 


We need not stop here. We can compute the tangent approximation to fat x, to 
obtain an even better approximation, x2, to the root of f (x) = 0, as in Figure 5, 


ys 


“5 


xo 
the tangent approximation 
at xy 


Figure 5. The tangent approximation at x,. 


The equation for x, can be deduced from Equation (2) as 
F(x) 
x2 =x, -—— 
mnie) 
Using this value of x, we can compute x3 using another tangent approximation at 


%2, and so on. Each term in the sequence xo, x,, x2, ... is computed using the 
recurrence relation 


I) 


fey (3) 


Xr+1 = Xp 


Example 2 
In Example 1 we looked for a root of the equation 
2x3 — 5x? — 5x +10=0 
near Xp = —1.5 and computed x; = — 1.478 7234. Continuing with this example, 


we build up the following table using Equation (3). (You may get slightly different 
results with your calculator.) 
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S(%,) = 2x3 — Sx,? — Sx, + 10 S'(&) = 6x? — 10x, — 5 


0 |-15 -0.5 23.5 — 1.478 7234 
1 | —1.4787234 —0.006 318 39 22.906 971 — 1.478 4476 
2 | —1.4784476 —0.000 001 06 22.899 319 — 1.478 4475 
3 | —1.4784475 —2x 107° 


After only three iterations f(x3) is so small that we can conclude that The precise value obtained for 
X3 = —1.478 4475 is a very good approximation to the root. J (x3) depends quite strongly 
e ee ¥ on which calculator is used 
There are two possible criteria for stopping the sequence: and how the calculation is 
(i) _Iff(x,) is very small, then x, is assumed to be very close to the root. Leegedah atlas 
(ii) If x,41 ~ X,, then Equation (3) gives 
LO) pias 5 L(%,) 
=x,- , which is equivalent to ——~ ~ 0. 
ae vt FG) 


Hence provided that f’(x,) is not too large we can conclude that x, is again a good 
approximation to the root. 

If we apply these criteria to Example 2, we find respectively: 

() f(s) = —2 x 10-8; 

(i) x; = —1.4784476, x3 = — 1.478 4475; 

and in either case we could conclude that there is no point in computing further 
terms in the sequence. 

We summarize the Newton-Raphson method in the following procedure box. 


The Newton-Raphson method 

Given an initial approximation xo close to the desired root of 
f(x) =0. 

Generate a sequence Xo, X;, X2, ... of successively better 
approximations using the recurrence relation 


_ f£) 
LY 


Xr = Xr 


This procedure can be terminated if 
either: (i) f (x,) becomes sufficiently small 


or: (ii) X41 = X,. 


To begin using the Newton-Raphson method we need an initial approximation Xo. 
If we do not have any information about x, the simplest thing to do is to sketch a 
graph of y = f(x), as we did for Example 1. The sketch need not be particularly 
accurate, and plotting half a dozen points will usually be sufficient to obtain crude 
approximations to the roots. In the following example, we suggest an alternative 
curve sketching method which sometimes can reduce the effort required to find a 
value for Xo still further. 


Example 3 
Find all the roots of 
x—2sinx=0 


using the Newton-Raphson method. Quote your answers to five significant figures. 
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Solution 


To use the Newton-Raphson method I need a value for xo. To obtain this I could 
sketch a graph of 


S (x) =x — 2sinx 


and look for the roots. However, it is easier in this case to draw graphs of y = x 
and y = 2 sin x and to look for points at which the two graphs intersect, since at 
these points we have x = 2 sin x. 


We can deduce that there can be no roots if |x| > 2, since —1 < sin x < 1, so that 
—2< 2sin x < 2 and the only possible roots of x = 2sin x are in the interval 
-2<x<2 


a 


Figure 6. Graphs of y = x and y = 2 sin x for -2<x <2 
From the graphs in-Figure 6 we can see that there are three roots. The first is at 
x = 0, while the other two roots are approximately —1.9 and 1.9. 
The Newton-Raphson recurrence relation is 
_ £65) 
L'%,) 
In this case, f(x) = x — 2 sin x and f’(x) = 1 — 2cos x. Hence the recurrence 
relation is 


Xe+1 = Xr 


x, — 2sinx, 


x, ——————————., 
rere" 1 = 2c0sx, 


With xo = 1.9 I obtained the following sequence, using my calculator. 


X= 19 
x, = 18955059 
X2 = 1.895 4943 (Uf (x2) = —9 x 10719) 


X3 = 18954943 (f(x3) = —4 x 10719) 


After only three iterations the values of x, agree to eight figures and f(x) is 
exceedingly small. 


Thus the root to five significant figures is 1.8955. 
Similarly the third root is at — 1.8955. Hence the three roots are 0 and +1.8955. 


In the above example the sequence converged very rapidly to the required root. 
The question we might ask is whether this will always be so. In the following piece 
of theory we will explain this rapid convergence and also gain some insight into 
how accurately we need to determine Xo, the initial guess at the root. 


Since we know that x = 0 is 
the exact value of a root, we 
do not need to use the 
Newton-Raphson method for 
this one. 


You may get slightly different 
results using your calculator. 


We know from basic 
trigonometry that 

sin (—x) = —sin x, so that if « 
is a root of x = 2 sinx then so 
is —a. 
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The recurrence relation for the Newton-Raphson method given by Equation (3) 
can be rearranged as 


O=S(%,) + (+1 — ¥)F (4) 


Let us compare Equation (4) with the first Taylor expansion about x, with 
remainder (page 11). Assuming that f(x) satisfies the conditions of Taylor's 
Theorem (with n = 1), we have 


I(&) =I Or) + % = He) fr) + HOE — HPI") (5) 
where c, lies between x, and x. 


If p is a root of the equation f(x) = 0, so that 


S(p) = 0, 
then we can substitute p for x in Equation (5) to get 
O =f (x) + (p — xf) + Ho — xP Feo) (6) 


where c, lies between x, and p. 
Subtracting Equation (6) from Equation (4) gives 
0 = (a1 — PIS) — Ho — xP Feo 


TA) 7 
Fe) ay 
Look at this last equation, (7), as this gives the secret of the success of the 
Newton-Raphson method. The left-hand side of the equation, x,,; — p, is just the 
error in the (r + 1)th approximation x,,, to the root p. Equation (7) relates this 
error to the error in x,. It shows that the error in x,,, is proportional to the 
square of the error in x,; so'if the error in x, is fairly small, the error in x,,, will 
be much smaller. As soon as we are reasonably close to the root, the errors 
diminish very rapidly. 


Xre1 — P= Hp — x)? 


In order to say more precisely when the method will converge, using Equation (7), 
we need to know some values for 


Ley) 

2f'(x,) 
which are generally not available. However, we can say that provided the 
conditions of Taylor’s Theorem are met and the equation f(x) = 0 does not have 
a root near p, then the sequence of iterates is guaranteed to converge for a suitable 
starting value xo. At the beginning of Subsection 2.1 I stated that we would 
assume that the function did not behave as in Figure 2. This figure shows a 
situation in which the equation f’(x) = 0 does have a root near p, so that the 
Newton-Raphson method would not work well. 


Exercise 1 
Find all three roots, to 5 significant figures, of the cubic equation 
x3 + 13.5x? + 40x + 16.67 =0 
using the Newton-Raphson method. Figure 7 shows a sketch of the graph. 
[Solution on p. 50} 


The results of this exercise will 
be used in Unit 24, Normal 
modes, 


If you are short of time, just 
calculate one of the roots. 
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gaat 3527 +40x +16-67 


—}— 
8 7 -6 -5 


-9 -4 


y(C-l0)=-35.5 


Figure 7. Sketch of y = x? + 13.5x? + 40x + 16.67. 


Exercise 2 
(i) Show that the roots of the equation 
xe*-0.1=0 
are at the intersections of the graphs of y = e* and y = 10x, and hence find rough 
approximations to them. 


(ii) Use the Newton-Raphson method to determine these roots to 6 decimal places. 
(Solution on p. 51) 


2.2 Taylor series methods for differential equations 
(Tape Subsection) 


The second application of Taylor polynomials that we look at is in the numerical 
solution of differential equations of the form 
y' =m(x, y) with y(xo) given You saw many differential 
uations of form in 


where m(x, y) is some given formula which may contain both x and y. init 2. 


For example, we may be asked to find an approximate solution to the differential 
equation 


Yaery 


for values of x such that 0 < x < 1, given that y = 0 when x = 0. This is the 
example we shall use in the tape frames which follow. 


(In this particular case we could compute an analytical solution using the The true solution is 
integrating factor method, but we use this simple example here to illustrate the J=e => 1s 
numerical method.) 


You have already met one method, Euler’s method, in Unit 2 for solving such 
problems numerically. The methods we shall develop here are called Taylor series 


methods. In general Taylor series methods will give much more accurate results 
than Euler’s method. 


Start the tape now. 
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Taylor polynomials 


Differential equation: = y' = x+y 
Initial condition + yeQ)=0 
PROBLEM : Use @ Ind Taylor polynomial approximation for 
Y about O to get an approximate value for y ak 2 =O! 
p(=) = y(0) + sey'(0) + Ya x*y"(0) (2) 
para: y(o) 


yo-L— 
vo] 


Computing the second derivative 
The differential equation (1) is 


y= ty 


"2 d(y)ed 
Als) a8 


Thus y= [ 
Hence yo) = [_____ 
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6 A new problem 
Differential equation : y zxtry (4) 
Initial_conditioa: = y (0:1)= 0-005 


PROBLEM: Use 8 2nd Taylor polynomial about O to 
get an approximate value for y a x= OD 


p(=) = y(Ol)+ (x-o1) y! (0:1) + Ya(x-0'1)*y"(01) (5) 
yo) LL] 


yo) LJ 
y"(o)= [7] 


ra] 


a 
sourion: pos)sL_______ | 
ed 


C4} The general problem 
Differential equation: y' = m (2,4) 


Initial condition : 4(=,) =o 


PROBLEM: Use @ 2nd Taylor polynomial about =<, to 
get an approximate value for 4 at c= Zee) 


P= y(2r) + (%- =r) Y!(%e) + Ya(mme,) y=) (6) 
where yy (x)= m (=..4(=)) (6a) 
and ye): go (mal) ) 


SOLUTION : p(z,,,)= 4(™) + (7 r) (=) + are =Jy',) (7) 
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Simplification 
Write y, for y(2r), uy, for y(=,) andy". for y"(x,) 
Write h = 2¢,,,-2, as the step size 


Equation (7) becomes 


P (ret) = Yet hy, + ¥ah*ye 
where yp =m (2,, 4-) 


and y = dm (=,,y,) 
doc 


Reworking the problem in Frame | 
Differential equation : y' ad | 
Initial condition : y(0) =O 

Taylor series method of onder 2 with h=O°l 
p(=,) = Yoethy, + Yah ye" (9) 
where ys = m(%,>4.) = + 4o* [] a) 
ard yf = (oes y.)= r=] 

a 


P(e!) 


Notation for the approximate solution 
Given Y,, we want values for y,,4.)44) °°, Yn 
Approximations : Given Y, = ye we actually compute NARA Ae A 
The Taylor series method of order 2 uses the recurrence relation 
to Eee +hY, + Vahty,” (10) 
where Ye = m (-<7,Ya) 


and" = dem (<2,,Y¢) 
da 
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A _tupicel problem 


Compute the approximate solution to the differential 
equation 


‘ 


Yr sry, 


with (0) =O, at <= Ol oa) 1 using the Taylor 
series méthod of order 2. 


ie ~ 
x, Y, | ied ea eect 


0 1 0.005 

0.005 0.021025 
0.021025 
0.049 23263 
0.09090205 | 0.49090205 | 1.49090205 0.147446 76 
0.147446 76 | 0.64744676 | 1.647446 76 0.220428 66 
0.220428 66 | 0.82042866 | 182042866 0.311573 66 
0.31157366 | 1.01157366 |2.01157366 0.422788 89 
0.42278889 | 1.22278889 | 2.22278889 0.55618171 
055618171 | 1.45618171 | 2.45618171 0.71408079 
0.71408079 
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Taylor Series method of order 3 


Compute the approximate solution bo the differential equation 


Giaaerys 
with 4(0) =O, at X=O1, 0-2,-+, 1, using the Taylor Series 
method of order 3. 


Cech 


warastas ) 
a 


¥, Y, =x,+7, 
7 

0 0 0.005 16667 
0.005 16667 | 0.105 16667 1.105 16667 0.02139335 


002139335") | ees 
0.049843 22 

0.09180173 | 0.49180173 | 1.49180173 0.148689 54 
0.148 68954 | 0.64868954 | 1.64868954 0.22207672 
0.22207672 | 0.82207672 | 1.82207672 0.313 69845 
0.31369845 | 1.01369845 | 2.01369845 0.4254724 
0.4254724 | 1.2254724 | 2.2254724 0,55951791 
0.55951791 | 1.45951791 | 2.45951791 0.71817721 
0.71817721 


Corr anewn 


Ss 
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We summarize the procedure for the Taylor series method of order 2 in the 
following procedure box. 


Taylor series method of order 2 
1. You are given a differential equation 


y! = m(x, y) 
and a value for y(xo). 


2. Differentiate the expression for m(x, y) to obtain or, y) 
remembering that y is a function of x. 


3. To apply the Taylor series method of order 2, choose a step 
size h, and calculate Y,, Y, ... from 


Yer=¥%+hy! + dh?y,” 
where 
¥, = m(x,, ¥,) 
dm 
i me ¥,) 
Yo = y(xo) and x, = Xo + rh. 


4. Y, is an approximation to y,. 


The procedure for using the Taylor series method of order n can be stated, in a 
similar way, as follows. 


Taylor series method of order n 
1. You are given a differential equation 


y’ =m(x, y) 
and a value for (xo). 


2. Differentiate the expression for m(x,y) n — 1 times to obtain 


F dm dm a'm 
expressions for Praaodd Be ys ase Bert Ow) 
remembering that y is a function of x. 


3. To apply the Taylor series method of order n, choose a step 
size h, and calculate Y;, Y2, ... from 
Pris enu ae 
Yi =%+hY +7 i Seclian serp 
where 


¥/ = m(x, ¥) 
» dm 
Y= FOr &) 


a 'm 
Pa 
Yo = y(Xo) and X, = Xo + rh. 


¥"= (&, %) 


4, Y, is an approximation to y,. 


There is a decision to be made here whether to use a higher-order Taylor series 
method with a larger step size or a lower-order Taylor series method with smaller 
step size in order to achieve a required accuracy. Reducing the step size can reduce 
the error significantly, but requires correspondingly more computational effort to 


The notation = (x, ¥,) means 


that we evaluate me.) at 


x=x,and y= ¥, 
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obtain approximations over a range of x-values. Adding an extra term in the 
Taylor series would probably not involve significantly more work provided that 
the derivative were easy to obtain. We shall say more about such comparisons in 
the next unit. 


One of the main drawbacks of the Taylor series methods is that they involve the 
use of higher derivatives which may be very complicated expressions. For example, 
you might like to convince yourself that the implementation of the Taylor series 
method of order 4 for the differential equation 


ya Very 
is by no means easy. (Do not spend more than a couple of minutes doing this.) 
In addition the work involved in inputting these higher derivatives into a 
computer or programmable calculator may be prohibitive. Consequently in the 
next unit we will look at other methods, which do not require these higher 
derivatives, 


Exercise 3 


Write down the recurrence relation for the Taylor series method of order 1. Does this 
method have another name? 


Exercise 4 

Use the Taylor series method of order 2 to solve the differential equation 
y =3y+sinx — with y(0) =0. 

Use h = 0.2 to compute approximate solutions at x = 0.2 and x = 0.4. 


Exercise 5 


Write down, the recurrence relation and expressions for ¥,’ (in terms of x, and ¥,), ¥," {in 
terms of x,, ¥, and ¥,) and ¥,” (in terms of x,, ¥,, ¥; and ¥,") required to use the Taylor 
series method of order 3 to solve the differential equation 


y =siny with y(0) = 1. 
Taking h = 0.1, compute the approximate solution at x = 0.1. 
[Solutions to Exercises 3-5 on p. 51) 


Summary of Section 2 


1. A root of the equation 
f(x) =0 
is any value of x for which the equation is satisfied. 


2. The Newton-Raphson method, given an initial approximation xo close to the 
root of 


L(x) =0, 
can be used to generate a sequence Xo, X;, X2, ... of successively better 
approximations to the root using the recurrence relation 
_ £0) 
f'%) 
The error in x,,, is proportional to the square of the error in x,, and so if the 
sequence converges it does so very rapidly. 


3. The Taylor series method of order 2 for solving first-order differential equations 
is given on page 27. 


Xr = Xp 


4. The Taylor series method of order n for solving first-order differential equations 
is given on page 27. 


5. The selection of the appropriate order of the Taylor series method and the 
choice of step size are difficult to determine for a given differential equation. The 
drawback of the method is that higher derivatives may be complicated expressions 
which are tedious to work out and time-consuming to input to a computer or 
programmable calculator. 
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3 The catenary (Television Section) 


To illustrate how some of the techniques of this unit are used in practice, we are rope 

going to look at a particular problem which arises out of the study of catenaries. ee al 
A catenary is the shape taken by a freely hanging rope or cable fixed at both ends, 

as shown in Figure 1. Figure 1. A catenary. 


Here is a statement of the problem we are going to tackle: Take a length L of 
Tope (where L is greater than 2m) and suspend it between two points 2m apart in 
the same horizontal plane. What is the displacement d of the middle of the rope 
below this plane? (L and d are shown in Figure 2.) 


Figure 2. Television problem is to find d given L. 


This problem has applications in the design of telephone cables or electricity 
cables, in that temperature changes cause the cables to increase or decrease in 
length. 


It would be useful if we could find a formula expressing the displacement d in 
terms of the length L, because it would then be possible to determine d by 
substituting the given value of L in that formula. Unfortunately such a formula is 
impossible to find. However, we shall see that the Newton-Raphson method can 
be used to determine d for any given value of L. It is also possible to use Taylor 
polynomial approximations to determine an approximate relationship between d 
and L which gives good answers for small displacements. To tackle this problem 
we need two pieces of information: 


(i) the equation expressing the catenary as a function of x; 
(ii) a formula for the length of a curve. 


These two topics are the subject of the next subsection, which you should read 
before watching the television programme. 
3.1 Hyperbolic functions 


The modelling of the catenary, which assumes that the mass per unit length is 
constant along the length of the rope, leads to a second-order differential equation. 


Unfortunately we do not have time to go into the details here. For a suitable They are given, for example, 
choice of the x- and y-axes the equation of the catenary can be written as in Differential Equations with 
Applications and Historical 
1 [e*+e"= Notes by G. F. Simmons, 
yes — (1) McGraw-Hill, 1972. 


where a is a constant whose value is chosen to select the particular catenary we 
want from a whole family of catenaries. (For each length L of rope we will have a 
different value of a,) The suitable choice of axes (see Figure 3) is such that 


y0)=0 
1 
and (0) ==. 
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y(0)=0 


a4 
-1 0 1 


Figure 3. Axes for the catenary in the case a = 1. 


The expression in brackets in Equation (1) appears so often that it is given a 
name. We define the hyperbolic cosine of the number ax as 


a+e™ 
— = 
The corresponding hyperbolic sine is defined as 


cosh ax = 


sinh ax = 


The equation of the catenary can now be written as 
y= cosh ax. (2) 


The hyperbolic functions have some very interesting properties, which we will use 
in the television programme. These properties are very similar to the properties of 
ordinary sines and cosines. In fact, to distinguish between the two types of 
function we call the well-known sines and cosines circular functions. This name 
arises because any point on a unit circle can be defined in terms of the cosine and 
sine of the angle @ shown in Figure 4. 


Since the equation of the unit circle is given by 

e+y=1, 
the well-known equation 

cos? @ + sin? = 1 
shows that the point (cos 0, sin @) lies on this circle. 
To see why we call cosh and sinh hyperbolic functions, consider the curve, called a 
hyperbola, whose equation is 

xt-yral. (3) 
I shall show that the point whose coordinates are (cosh @, sinh @) lies on the 


hyperbola, i.e. that it satisfies Equation (3). The number can take any value, but 
its geometical significance is less obvious than for the circle. 


Now 


“ — 
cosh 0+ sinh 9 = [°* J+ < )-¢ 


2 
and 
f lP+e* l?—e® = 
cosh 6 — sinh 6 = | - ae = 
Thus 


cosh?@ — sinh?@ = (cosh @ + sinh @) (cosh @ — sinh 0) 
=ee*=1 


as required. 


The symbol cosh is 
pronounced as ‘cosh’. 


The symbol sinh is 
pronounced as ‘shine’. 


{cos @,sin0) 


Figure 4. Any point on the 
unit circle can be expressed in 
terms of 8. 


(cosh0, sinh0) 


Figure 5. Any point on the 
hyperbola x? — y? = 1 can be 
expressed in terms of 0. 
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Note the similarity between the result 
cosh?@ — sinh?0 = 1 

for hyperbolic functions and the result 
cos?0 + sin?@ = 1 

for circular functions. 


The following table lists some properties of hyperbolic functions together with the 
corresponding properties of circular functions. They can be derived using the 
definitions of cosh x and sinh x in terms of exponential functions. If you have time 
you might verify some of these properties. You are not expected to memorize 
these; they will be included in the Handbook. 


Hyperbolic functions Circular functions 
e+e* e* +e7% 
cosh x = os = 
A ee” : 
sinh x = sin x = 


cosh?x — sinh?x = 1 cos?x + sin?x = 1 


d : d 5 
a (Cosh x) = sinh x a (008 x) = —sin x 


are oe 
ae inh x) = cosh x aon x) = cos x 


ext xf ef 
cosh cele rhea hs conse a 
‘ah SoCo ee ioe eee 
sinl Sdikalr 25-3 / at be ahsielsSoer ns 


NO OO ——————————————E 
The last two properties in the table can also be derived using Taylor’s Theorem, 
given in Subsection 1.3. For example, if f(x) = cosh x, the conditions of Taylor's 
Theorem are satisfied for any interval [a, b] and we may write the 2nth Taylor 
expansion about 0, with remainder, as 


2 
cosh x = cosh 0+ x sinh 0 + 5 cosh 0 + «++ 


x? aati ; 
+ Giro 0+ @an cy (4) 


where c, lies between 0 and x. 


Now cosh 0 = 1 and sinh 0 = 0. Hence Equation (4) simplifies to 


2n xaetl 


Bip an es 
(2n)! (Qn+1)! 


As n becomes large we obtain the required expression for cosh x, since for any x 
the remainder term approaches zero for large n. (We will not prove this, as it 
would take us too far from the main subject of the unit.) 


2 4 
cosh x =14+ 545 +- sinh ¢,. 


2! 


In Exercise 3 you will be asked to obtain a formula for the length of a catenary. 
To do this you will need to use the following result. 


Length of a curve 
The length L of the curve y = f(x) between xo and x, is given by 


‘These first two expressions for 
the circular functions were 
given in Unit 5, Subsection 
42. 


You can verify this result in 
the case of the straight line 
y=ax 
between xp = 0 and x, = 1 by 
showing that 
L=/1+a. 
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For most curves a numerical method would have to be used to evaluate L, but in 
the case of a catenary we can determine L exactly. 


Exercise 1 
Show that: 
(i) cosh (—x) = cosh x 
sinh (—x) = —sinh x 
(ii) 4 cosh x) =sinh x 
A ion x) = cosh x. 
‘Compare these results with the corresponding properties of sines and cosines. 


Exercise 2 
Use Taylor’s Theorem to derive the Taylor series for sinh x as 


. xe ae 
ex+ytatate. 
sinh x x+5 TI 


Exercise 3 
Determine the length of the catenary, between x = —1 and x = 1, given by 
1 
y= jo08h ax. 


[Solutions to Exercises 1-3 on pp. 51-2] 


Now watch the television programme: *Catenaries — numerical approximation’. 


3.2. The rope problem 


As I stated in the introduction to Section 3, the problem considered in the 
television programme is to find the displacement d at the middle of a length L of 
rope suspended between points two metres apart as in Figure 2. 


In Subsection 3.1 I stated that the equation of a catenary is given by Equation (2) 
as 


y= cosh ax. (5) 


y=2 cosh ax 


+ ae 
=! 0 1* 
Figure 6. The catenary. 
To obtain an expression for d we note that the height of the catenary above the x- 
axis at x= 1 is d-+ + (ee Figure 6) ic. 
i 
yl) =d+-. 
a 
Using the equation for the catenary (5) we have 


Hl) eet am da 
a a 
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Hence 

d= + (cosh a1), (6) 
This gives d in terms of a, whereas we really require an expression for d in terms 
of L. 


From Exercise 3 we know that the length of the catenary is given by 
L=?2sinh a 7 
a 


Equations (6) and (7) define d and L in terms of a, but unfortunately there is no 
way to establish a direct formula for d in terms of L. It is possible to use these 
formulae, with the help of the Newton-Raphson method, although it does involve 
a lot of work. 

Example 1 

Determine the displacement d if L = 2.2m. 

Solution 

From Equation (7) we have 


22 aiohi a (8) 
a 


Rearranging this equation gives a as a root of the non-linear equation 
f(a) =sinha—1.1a=0. (9) 


From the sketch graph of f (a) in Figure 7 we can see that there are two non- 
negative roots for Equation (9). However, the root a = 0 is not a solution of 
Equation (8) and the value of a we require is approximately a = 0.76. 


Figure 7, Sketch of f(a) = sinh a — L.la. 


The Newton-Raphson method, derived in Subsection 2.1, determines a root of the 
equation f(x) = 0 using the recurrence relation 


with Xo given, 


For our function of a we have 
f(a) =sinha-1.la and = f(a) = cosh a— L.1. 
Hence the recurrence relation is 


sinh a, — 1.1 +) 


Dey Oy — | cosh a, — 1.1 | 
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With ao = 0.76 we have 


a, = 0.763 424 64 
a = 0.763 40080 
= 3. 


After only 3 iterations the sequence has converged to the required root. With 
a= 0.763400 80 we can now use Equation (6) to evaluate d as 


d= *(cosh a — 1) = 0.400 601 6. 


Thus when L = 2.2m, the displacement at the centre of the rope is given as 
d =0.401 m correct to 3 decimal places. 


In the above example we have managed to find d for a given length L. However, 
for each value of L we have to determine an initial value ao, use the Newton- 
Raphson method, and then substitute the computed value of a into Equation (6) 
to determine the displacement. 


In the second part of the television programme we look for a simple approximate 
relationship between d and L using Taylor polynomials, which works well for 
small displacements d. The advantage of a simple relationship is that it will be 
easier to make qualitative predictions about the behaviour of the rope for small 
displacements. Thus, although it is possible to obtain a more accurate 
approximate relationship between d and L, the simple one may be more useful in 
practice. 
From the properties of hyperbolic functions we know that 

4 


(ae tease 
cosha=l+o +a tate: 


Hence the 2nd Taylor polynomial about 0 for cosh a is 


cosh a = 1 + 4a. Because the a? term is zero, 
e . Be ers 7 ae =: this is actually the 3rd Taylor 
While a is reasonably close to zero this will provide a good approximation to polynomial as well as the 2nd. 


cosh a. In Equation (6) we have 
1 
d= q(cosh a-1) 


= 1 + 4a? —1) (using the 2nd (or 3rd) Taylor 
ks polynomial about 0 for cosh a) 
=ha. (10) 


Similarly we have 
3 a 


inh a 
sinha=a+ta p+ opto" 


and the 3rd Taylor polynomial about 0 for sinh a is 
sinh a~a + a’. 


Thus the equation for L is 
L= 2 inh a 
a 
2 3 
=7G + $a’) 


=2+4a?. (11) 


Now Equation (11) can be rearranged to give a in terms of Las 


ax ./3(L—2). 


Hence using Equation (10) the approximate displacement d* is given by 


d* = },/3(L — 2) = \/4(3L — 6). (12) 
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Now we have an explicit expression for d* for any value of L. For example, if 
L=2.2 we have 


d* = ,/0.15 = 0.387, 
which is quite close to the more accurate value d ~ 0.401 m, obtained in Example 1. 


Comparing the approximate results, for different values of L, with the true results 
in Figure 8, we can see that for values of L slightly greater than 2 the 
approximation is very good, but as L increases the approximation gets worse. 


d 
v true values. 
oes . 
0S approximate solution 
0 eer 


== = =. se 
25 2.6 27 L 


Figure 8. Graph comparing the approximate displacements with the true displacements, 


— aly 
2.0 21 22 23 24 


We conclude that our approximate method is a very useful way of determining the 
approximate displacement d* when d is small. Furthermore, the approximate 
relationship tells us something we might have otherwise missed. For L slightly 
greater than 2 there is a very sharp increase in the displacement. For example, if 
L=2.001 we find that the displacement d ~ 0.027, ie. if L is increased by 1mm 
from the zero-displacement position then the displacement will be 2.7cm, 27 times 
as big. Hence the problem of finding d given L is absolutely ill-conditioned for L 
close to 2m. The ill-conditioning is demonstrated at the end of the television 
programme, illustrating that approximate relationships can give us information 
that would otherwise be difficult to deduce. 


Exercise 4 

Use the Newton-Raphson method to determine the displacement d, to an accuracy of 6 
significant figures, when L = 2.08 m. 

Exercise 5 

Express a in terms of L if the Sth Taylor polynomial about a = 0 is used to approximate 
sinh a in Equation (7). 

Exercise 6 


Use the approximation derived in Exercise 5 together with Equation (6) to find the 
approximate displacement if L = 2.08m. 


[Solutions to Exercises 4-6 on pp. 52-3] 


Summary of Section 3 


1. The equation modelling a freely hanging rope suspended between two points is 


ve) a al 
a 2 


where a is a parameter which depends on (among other things) the length of the 
rope, and cosh is the hyperbolic cosine function. The curve with this equation is 
called a catenary. 


Each of the true values was 
computed using Equations (6) 
and (7) and the Newton- 
Raphson method. 
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and 


have many properties which are similar to those of the circular functions sin x and 
cos x. A list of these properties is given on page 31. 


3. The length of a curve between x» and x, is given by 


1: fdy\? 
L= is 1+ (?) dx. 
LN, ax 
In most cases a numerical method would be required to evaluate this integral. 
4. The displacement d for a freely hanging rope of length L suspended between 
two points 2m apart is given by 
dm (Cosh a —1) 
a 
where 
L= 2 inh a. 
a 


For a given value of L, a can be determined using the Newton-Raphson method 
and this value can be used to determine d. An approximate relationship between d 
and L, obtained using the 3rd Taylor polynomials about a = 0 for cosha and 
sinh a, is given by 


d* = ,/}(3L — 6). 


For 2 < L < 2.4 this gives approximations within 7% of the true solution. 


4 Interpolation and integration 


4.1 Interpolation 


In Section 1 we derived the Taylor polynomials of degree n by assuming that we 
knew or could determine n derivatives of a function at a particular point a. 
However, in many problems the values of the derivatives are not known, but the 
data is given at a number of points Xo, X1, X2, X3, ... , X». Diagrammatically we 
could represent this as in Figure 1. 


S(x9) 
x s(x) 
x) fx) 
x 
T T = 
a a ae Xn x 


Figure 1. Data given at n + 1 points xo, x1, -.- Xe 


It is possible to show that there is just one polynomial of degree n which takes the 
desired value of f at each of the n + 1 points xo, X1, -.., Xx 
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Once we have determined this polynomial we can use it to approximate the 
function at other values of x. This process of determining the approximate value of 
the function at a point x between xo and x, is called interpolation. (The process of 
using the polynomial to determine approximate values outside this range is called 
extrapolation. Extrapolation is much less accurate than interpolation and we will 
not consider this topic here.) One method of computing the coefficients of the 
interpolating polynomial of degree n, p(x), given as 

P(X) = Qo + ayx +... + a,x", 


is to set up a system of equations. Since p(x;) = f(x,) for i= 0, 1, ..., n, we can 
write 

P(Xo) = ao + ay Xo +... + aX" =f (Xo) 

P(t) = do + 1X1 +... + ays" =f) 


P(Xp) = Ao + GyXq +... + GpXq" = f (Xp) 
This is just a system of linear equations, like those you met in Unit 9, Simultaneous 


linear algebraic equations, for the n + 1 unknowns do, dj, ..., d,. In matrix form 
the problem is to find ao, a), ..., a, which satisfy the equations 


A) xq. oxor Xo" | [ao I (Xo) 

lox x"! fal | fe) 

LXe Xe xn" J Lay, () 
Example 1 


Determine the interpolating quadratic polynomial which passes through the three 
points (0, 1), (1, 2) and (2, 4). 


Solution 
Let the interpolating polynomial be given by 
P(X) = ao + a,x + apx?. 
The three equations corresponding to the three data points are 
(0) = do + a, (0) + a2(0*) = 1 
P(1) = do + ay(1) + a (1?) =2 
P(2) = do + a;(2) + a2(27) = 4 


giving 
ay =1 
ago+ a+ a,=2 
ao + 2a, + 4a, =4 
ie. 


12 4}/a,| |4 


Using Gaussian elimination, as described in Unit 9, we reduce this matrix to its 
upper triangular form as 


1 0 Ojfa} {1 
0 1 tifa, j=]1 
0 0 2})a,} {1 


whence 


Note that in this problem we 
know the values of xo, x, ... 
%y and f (Xo) f (1)... Xm 
and we wish to determine 
9, y+ Aye 
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Thus the interpolating polynomial is 

p(x) = 1+ 4x + 4x?. 
Given a large set of points it is not always desirable to compute the interpolating 
polynomial of corresponding degree. This process can be very time-consuming, 
prone to the build-up of rounding errors in computation, and the resulting 
polynomial tedious to evaluate using a calculator. More frequently, given a value 
of x at which we need to approximate the function, we construct a polynomial 
using just a few of the data points. 


Example 2 
Here is a set of data obtained from four-figure tables for log, x. 


~ logsr%e 


Use an interpolating polynomial of degree 2 to obtain an approximation to 
log, (1.134). 


Solution 


For a quadratic polynomial we need to specify 3 points xo, x, and xz. We choose 
the three points closest to x = 1.134, namely 


Xo = 1.12, x, = 1.13, X2 = 1.14. 
If p(x) = ao + a,x + ax, we obtain the following equations for do, a; and a3: 
p(1.12) = ao + 1.12a, + (1.12)? a2 = 0.1133 
p(1.13) = do + 1.134, + (1.13)? a, = 0.1222 
p(1.14) = ay + 1.14a, + (1.14)? az = 0.1310 


1 1.12 1.2544 |}a9 0.1133 
1 1.13 1.2769 |} a, | =| 0.1222]. 
1 1.14 1.2996 || a2 0.1310 
Using the computer package SIMLIN,I obtained the solution 
ay = —1,51635, a, = 2.01508, a, = —0,50004. 
Hence 
p(x) = —1.51635 + 2.01508 x — 0.50004 x?, 
giving 
p(1.134) = 0.1257 to 4 significant figures. 


It is worth pointing out that the problem of evaluating the coefficients of p in the 
above example is highly ill-conditioned. To illustrate this I repeated the 
calculations with the right-hand side values correct to five significant figures. The 
new set of equations to solve is 


1 112 1.2544][ag] [0.11333 
1 143 1.2769] a,| =| 0.12222 |. 
1 114 1.2996||a;| | 0.13103 
SIMLIN this time gives 
dy = —1.38863, a, = 1.78908 and a, = — 0.40004, 
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so that 

P(x) = —1.388 63 + 1.78908 x — 0.40004 x?. 
This is completely different from the polynomial obtained in Example 2, and the 
difference is caused by very small changes in the data. The reason for this ill- 
conditioning is partly due to the fact that the rows of the matrix are almost 
linearly dependent (you can verify this if you wish!). However, if we use the new 
polynomial to evaluate an approximation to log, (1.134) we get 


p(1.134) = 0.12575 to 5 significant figures, 


and this answer is almost identical to the answer obtained in Example 2. The true 
solution (to 8 decimal places) is 


log. (1.134) = 0.125751 21, 


indicating that the answers we obtained are remarkably good. Although the 
polynomials are completely different, the difference in their values for 

1.12 < x < 1.14 is less than 0.00004, ic. approximately the same as the differences 
in the data! 


What has happened is that although the problem of computing the coefficients of 
the interpolating polynomial is ill-conditioned, the problem of determining the 
approximate value for log, (1.134) is well-conditioned. This phenomenon occurs 
very frequently in interpolation problems where the x-values are close together. 
You need to be aware that this can occur, although it almost always makes no 
difference to the final answer. 


Exercise 1 


Determine the quadratic whose graph passes through the three points (—1,2), (0,3) and 
(1,0). 


Exercise 2 
Given the data in Example 2, use an interpolating polynomial of degree 1 to compute an 
approximation to log, (1.134). 


(Solutions to Exercises 1 and 2 on p. 53) 


In Section 1, when we discussed Taylor polynomials, we stated Taylor's Theorem 
which could be used to determine the error in the approximation. In the same way 
there is a very useful theorem which can be used to determine the error in the 
interpolating polynomial. In this theorem we assume that all arithmetic has been 
carried out exactly and there are no rounding errors in the values 


SL (Xo) f 1-2 f On). 


Theorem 4.1 
You are given a function f(x) and an interpolating polynomial 
p(x) of degree < n such that 


P(x) =f(x)  fori=0,1,..., 0, 
where the numbers Xo, X;, ..., X, ate arranged in increasing order. 
Then the error, e(x) = p(x) — f(x), where x lies between xo and x,, 
is given by 
*O(c.) 

(n+ 1)! 


where c, is some number satisfying xo < c, < x» 


@(x) = —(% — x0)(x — x1) +++ (& — x, 


(ly 


You can see that the presence of the product (x — x9) (x — x)... (x — x,) ensures that 
e(x;) = 0 for i=0, 1,...,” 
as required, while the term 


£* P(e) 
(n+ 1)! 


is very similar to the corresponding term in Taylor's Theorem. 
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The following procedure uses Theorem 1 to compute an error bound at a particular 
value of x between xo and x,. 


Error bounds for interpolating polynomials 
1. Let p(x) be the polynomial of degree <n which interpolates the 
function f(x) at Xo, X1, ..., Xm Where x9 <X, <...<Xp 


2. The error e(x) = p(x) — f(x) is given by 
ax) = — (x — x0)(x — x1)... — io 
where c, lies between xy and x,. 
3. Determine a value for M such that 
If < M 
for all values of c between xo and x,. 
4. A bound on the error e(x) is then 


I = x0)(x — X1)..-(% — Xn)|M 
(n+ 1)! 


le(x)| < 


To see how this procedure can be used, we calculate an error bound in the 
following example. 

Example 3 

(0) Show that |1/c?| < 1 for all values of c between 1 and 2. 


(ii) The linear approximation to log, x whose graph passes through the points 
(1, 0) and (2, 0.6931) is given by 


P(x) = 0.6931x — 0.6931. The only error in the 

P c 2 calculation of p(x) is in giving 
Determine an error bound at x = 1.5 and compare this bound with the log, 2 to only 4 decimal 
actual error at x = 1.5. places. 


Solution 
(i) Since c > 1, we have 0 < I/c < 1, and so |1/c*| < 1. 


(ii) Using the procedure for computing error bounds, we note that with 
Xo = 1 and x, = 2 the error function is given by 


dfx) =~ — 1-2 


where c, lies between | and 2. 
Now f(x) = log, x, and so 
Sx) =1/x and f"(x) = —1/x?. 


Hence from part (i) an upper bound for |f”(c)| = |1/c?| where 
1 <c <2 is given by 


M=1. 
Thus |e(x)| < fen Ne 2H and at x = 1.5 we have 
pa.sy HS ES—2 gras 


Thus an error bound for this approximation is given by 0.125. The actual 
error is given by 


p(1.5) — log, 1.5 = 0,346 55 — 0.40547 
= 0.05892 


so that |e(1.5)| = 0.058 92. This is certainly less than the error bound 0.125. 
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Exercise 3 

@ Show that, for —1 <c < 1, we have |e] <e = 2.7183... 

(ii) The function f(x) = e* is to be approximated by the interpolating quadratic 
polynomial which interpolates f at x» = xX, =O and x, = 1. Determine a 
bound for the error in this quadratic at x = 0.5, ignoring rounding errors. (Note 
that you are not asked to compute the interpolating quadratic polynomial.) 

[Solution on p. 53] 


4.2 Methods of numerical integration 


In the last subsection we looked at methods of interpolation. In this final 
subsection we examine how these interpolation methods can be used to derive 
approximations to the integral 


t= [yoyex 


where f is some given function. 


By integrating the interpolating polynomial, which is a fairly straightforward 
Process, we can obtain formulae for the approximation of the integral. Such 
formulae are known as Newton-Cotes formulae. The methods derived here, Euler’s 
integration method, the trapezoidal method and Simpson’s method, will be used in 
the next unit to solve differential equations, but they are clearly useful in their own 
right as well. 


Euler’s method The reason for calling it 


A - Euler's method is that it i: 
Suppose we want to approximate the integral a aay rated ie Euler's " 


1 method for differential 
f S (x)dx. equations. See Unit 19, 
x0 
Figure 2 shows one simple approximation: we replace the function f(x) by the is 
constant f (xo) for all values of x between xo and x;. This is equivalent to using a : 
very crude ‘approximating polynomial’, given by function y'= f(r) 
P(x) =f (Xo). 
The integral 
1 
t= ['seodx 
xo 
is the area under the curve y = f(x) between xo and x,. We approximate this area *o "a 
by the area under the approximating polynomial as Figure 2. Euler's method. 


Tf neowte= [" focodde = (x1 — olf) 
: ; =hf (xo) where h = x; — Xo, 


ive, the area under the curve y = f(x) is approximated by the area of the rectangle 
of height f (xo) and width h. 


Clearly this is a very crude approximation, and Euler's method would rarely be 
used in this form. However, we can get more accurate results by using a composite 
form of the method. Suppose we wish to approximate the integral of f over the 
interval [a, b]. Then we divide this interval into n subintervals of equal width h 
and use Euler's method on each subinterval (see Figure 3). 


The approximation to the integral is then 


f Slxjdx = [ “Ple)de $+ ib Seside 


= hf (Xo) + hf (x1) + ++ + Af (Xp—1) 
= (Sf (Xo) + fea) + + F%q—-1))- (2) 
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\ 


Y, 


| 


Figure 3. Euler's method in composite form. 


Ye 


Z 


\ 


AX YY) 


Example 4 
Determine the approximate value of 


2 
I= f log, x dx 
1 


using the composite Euler’s method with h = 0.1. 
Solution 
Using Equation (2) with h = 0.1, we have 


I~ 0.1 (log.1 + log,1.1 + log,1.2 + -+- + log.1.9) 
= 0.35122. 


For this example the true solution is given by 


2 
f log, x dx = [x log.x — x]? 
1 


= 2log.2—-2+1 

= 0.386 29. 
The approximation is not particularly good, but from Figure 3 we can deduce that 
if we use a very small interval width h we can achieve much better accuracy. 


Eauler’s method in composite form 
1, You are asked to obtain a value of the integral 


I= [reads 


2. Choose values for h and n such that 
b=a+nh. 


3. Approximate J using Euler's method in composite form as 
T= Af (Xo) +f (1) + °° +f n-1)) 
where x, = a + rh, 


Exercise 4 


Use Euler’s method in composite form with h = 0.2 to approximate the integral of e* 
between x = 0 and x = 1. 


[Solution on p. 53] 

The trapezoidal method 

Consider the linear interpolating polynomial approximation to the function 
y =f (x) which passes through the data points (xo, f (xo)) and (x, f(x1)). 
The integral 


I= [oreo 
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is approximated by the area under the graph of the interpolating polynomial. For 
linear interpolation this area is a trapezium—a quadrilateral with two opposite 
sides parallel. The area of this trapezium is the product of the width h = x, — x9 
times the average height, which is }(f (xo) +f (x,)). The approximation to the 
integral J is thus 


1 = th( (x0) + f(x1)). (3) 


function » = f(x) 


A~ interpolating 
polynomial 
y= p(x) 


x x 


xo 
Figure 4, The trapezoidal method. 


This method of approximating the integral is known as the trapezoidal method. We 
would expect the trapezoidal method to give more accurate results than Euler’s 
method, since the trapezium fits the area under the curve better than the rectangle 
does. 


Example 5 
Determine an approximate value of 


2 
1=f log, x dx 
1 


using the trapezoidal method. 
Solution 
Equation (3) gives the approximation as 
I = Slog, 1 + log, 2) = 0.34657. 
The correct solution is given in Example 4 as J = 0.386 29. Even with this fairly 
crude method we have obtained a reasonable approximation to the integral. 


The trapezoidal method is most commonly used in composite form. That is, if we 
wish to approximate the integrand f over the interval [a,b], then we divide this 
interval into n subintervals of equal width and use the trapezoidal method on each 
subinterval, as shown in Figure 6. 


function y = f(x) 


\\ 


ako x) ae = x,=6 x 


Figure 6. The composite trapezoidal method. 


I= [revue = HCL (0) + F (xa) + HAF (x1) +f (2) + 
+ HAUS p—1) + Fn) 
= h(BS (Xo) + f(x1) +f (e2) + + fe) + 46%). (4) 


Figure 5. The trapezoidal 
method for Example 5. 
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Example 6 
Determine the approximate value of 
2 
i= f log, x dx 
1 
using the composite trapezoidal method with h = 0.25. 
Solution 
Equation (4) with f(x) = log-x and h = 0.25 gives 
I =~ 0.25 (Slog, 1 + log, 1.25 + log. 1.5 + log, 1.75 + }log, 2 Figure 7. The composite 
Glog. al me Hog. 2) trapezoidal method applied in 
= 0.38370. Example 6, 


This is almost correct to 2 decimal places, since the true solution is J = 0.386 29. 


‘The trapezoidal method in composite form 
1. You are asked to find a value of 


I= f f(x)dx. 


Choose values for h and n such that 
b=a+nh. 


Approximate / using the trapezoidal method in composite 
form as 


T= W(bf (Xo) +f %1) +f (%2) + +L %n-1) +4) (5) 


where x, =a+rh. 


Exercise 5 

Use the composite trapezoidal method (Equation (5)) with h = 0.2 to estimate the integral 
of e* between x = 0 and x = 1. Compare your answer with the true solution and the 
solution obtained in Exercise 4. 

(Solution on p. 53] 

Simpson’s method 

Our third method of numerical integration is known as Simpson’s method. It is 
derived by interpolating three equally spaced data points (xo, f(xo)), (x1, f(%1)) 
and (x2, f(x2)) by a quadratic polynomial and then integrating this quadratic, as 
in Figure 8. Since it can allow for the curvature of the graph of the function f, we 
could expect Simpson’s method to give much more accurate results than Euler's 
method or the trapezoidal method. 


function y = f(x), 


interpolati 
quadratic palynomial 


Figure 8. Simpson's method. 
To simplify the derivation of the formula for Simpson's method we make a change 
of variable, writing 


S=x-X, (6) 
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and we denote the quadratic interpolating polynomial, written in terms of the 
variable s, by 


q(s) = a + bs + cs? 


From Equation (6) and the conditions determining the interpolating polynomial 
we see that: 


when x = Xo we have s = —h and q(—h) =f (xo), 
when x = x, we have s = 0 and q(0) =/(x,), 
when x = x2 we have s = h and q(h) =f (x2). 


Hence the conditions determining the coefficients a, b, c in the polynomial q(s) are 


q(—h) =a — bh + ch? =f (xo) (7) 
q0)=a =f (x1) (8) 
q(h) =a + bh + ch? = f(x,). (9) 


Equation (8) gives a = f(x,). Thus Equations (7) and (9) can be rewritten as 
—bh + ch? =f (xo) —f (x1) 
bh + ch? =f (x2) —f(x1). 
These two equations can now be solved for b and c, giving 


p —f 62) =F (xo) 
2h 
and 


I (X2) = 2fl%1) +f (X0) 
2 P 


as) =f (x1) + 


(Aer feed [Fetes feo) re 


To approximate the integral of f(x) between xo and x, we write 


i Sloxydx = [neste 


= J atows 


2 [fro (ef). 


¢ (fe - es +f) ‘4 a 


= 2 
= [san [POS 


. (eva 2f Se) sy} 
2h? Silay 


= 2hflos) +0 +5 (Foex) — 2F 64) +F¢00)). 


Hence 
fF foade =F (ra) + 4F 001) + F60)) (10) 


Now that we have the formula for Simpson’s method, given by Equation (10), we 
can simply use this to approximate integrals. 


Do not spend a lot of time 
working through the 
derivation of Equation (10). 


When x = xo, s = —h, and 
when x = x2, s = h, Equation 
(6) gives ds = dx. 
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Example 7 
Determine an approximation to the integral 


2 
I= f log. x dx 
1 


using Simpson’s method. 
Solution 
To use Simpson's method we need two equal intervals of width h = }. From 
Equation (10) we have, since h = $(x2 — Xo) = 4, 
1 = } (logel + 4log.1.5 + log,2) 
= 0.38583 
and this is correct to three decimal places. 


Like Euler's method and the trapezoidal method, Simpson's method can be used 
in composite form. We divide the interval [a,b] into an even number of 
subintervals and approximate the integral as the area under the quadratic for each 
pair of subintervals. 


function 
=! yas) 


aux) %) x2 xy Xe-1 Xq_ 7b 
ee 2) + —+| 1 


Figure 9. The composite Simpson's method. 


[pecrte = 8 reve) +49 (es) + F00)) +5 (Fe) + AF 9) +60) + 
h 
+ 3 Gra) + 4S (%n-1) +f %n)) 


h 
= 3/0) + 4f (x1) + 2X2) + AF (Xs) + ZF (Xa) + °° 


+ 4f %n-1) + fn). ay 


The composite Simpson’s method is one of the most popular methods for 
numerical integration. We summarize the procedure in the following box. 


Simpson’s method in composite form 
1. You are asked to obtain a value for the integral 


I= [ree 


2. Choose values for h and n such that n is even and 
b=a+nh. 


3. Approximate J using Simpson's method in composite form as 


1 (feo) + 4fl0) + 2f 62) + Af es) + 2f ls) + 
+ 4f op=1) + (0) 


where x, = a+ rh. 


We have not drawn a figure in 
this case, as there would be no 
visible difference between the 
graph of y = log, x and the 
interpolating quadratic. 
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Example 8 
Determine an approximate value for the integral 


2 
1 log.x dx 
\ 


using the composite Simpson's method with h = 4. 
Solution 
With h = 4, Equation (11) gives 


T= 7; (loge! + 4log.1.25 + 2log,1.5 + 4log,1.75 + log.2) 
= 0,386 260 


and this solution is correct to four decimal places, 


Exercise 6 
(i) Use Simpson's method to obtain an approximation to the integral of e* between x = 0 
and x = 1. 


(ti) Use the composite Simpson’s method for the same integral with h = 4 and h = §. 
Check the accuracy of your solutions. 


Exercise 7 
The length of the curve y = sin x between 0 and ; is given by 


(2 
L={ (1 + cos? x)!/? dx, See Subsection 3.1. 
0 


Use the composite Simpson’s method with h = ; to find the approximate length. 


[Solutions to Exercises 6 and 7 on p. 53) 


Summary of Section 4 


1. Given a set of data (xo, f(xo)), (X1.fl%1)), --- 5 (Sm S%m))s interpolation refers to 
any method of estimating the value of f(x) at an intermediate point x by 
constructing an interpolating function whose graph passes through these 
points and then evaluating this function at x. 


2. To construct an interpolating polynomial p of degree <n given by 
P(X) = do + ayxX + ayx? + +--+ a,x" 


which passes through the points (xo, (xo), (x1, f(%1)),. + Xn SO), we solve 
the n + 1 linear equations for do,...,a, given by 


P(x;) = f(x), 1=0,1,....m. 


3. The interpolation error for a polynomial of degree n is given in Theorem 4.1. A 
bound for this error can be found using the procedure on page 40. 


4. Euler's method for approximating the integral of f between xq and x, is given 
by 


[fords = hy xo) 


Euler’s method in composite form for approximating the integral of f between 
Xo and x, is given by 


[pera = b Fo) +f) + Fe) +--+ F054) 


whereh=x;.,-x,  i=0,...,n—1. 


6. 


5 
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The trapezoidal method for approximating the integral of f between xo and x, 
is given by 
- A 
fF Feadx =F eo) +f) 
xo 


where h = x; — Xo. 


The composite trapezoidal method for approximating the integral of f between 
Xo and x, is given by 


J feoads = hf 00) +S) + Fbe0) + + Fe 0) + AF 0) 
where h = x41 — Xi i=Q....n-L 


‘Simpson’s method for approximating the integral of f between xo and x; is 
given by 


J soot ~ E00) + 4f 0) + F003) 


where h = x2 — x, =X1 — Xo. 
The composite Simpson’s method for approximating the integral of f between xo 
and x, (n even) is given by 
J ses = 5 G0) + Af 1) + 24 2) + AF 9) + 20+ 27 0) 
+4 Sl%q-1) +S (%0)) 
where h = x44 — Xi, i=0,...,.0—1. 


Simpson's method is usually much more accurate than the trapezoidal method or 
Euler's method. 


End of unit exercises 


Exercise 1 
(i) | Show that the maximum value of 2 tan csec*c is 0.20 to 2 decimal places for Hint: 
O<c<Ol. sec?x = 1 + tan?x 
(ii) Determine the linear Taylor polynomial about 0 for tan x and give an error bound Sig = sec?x 
‘on the approximation of tan x by this polynomial for 0 < x < 0.1. 
Exercise 2 sce x = sec x tan x, 


Determine the two positive roots of the equation 


f(x) =e -3x=0 


to 5 decimal places. 


Exercise 3 


Use the Taylor series method of order 2 with h = 0.1 to obtain approximations to the 
solution of the differential equation 


y=xy+1 with yQ)=1 


at x = 0.1 and x = 0.2. Give your approximations to 3 decimal places. 
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Exercise 4 
(i) You are given just the following information about a function f: 
S(1) = 1 and f (24) = 14. 
Use an interpolating polynomial of degree one to approximate the function at 
44, 


x= 144. 
(ii) Verify that for 1 < ¢ < 2} we have [fe73?| <4. 
(iii) The function f in this problem is actually given by 
Sx) =./x. 
Determine an error bound for the approximation at x = 1.44 and compare this 
bound with the actual error. 
Exercise 5 
Use Simpson's method in composite form with h = 0.25 to approximate the integral of 
tan x between 0 and 1. 
Exercise 6 


Use Simpson's method in composite form with h = 0.25 to determine the approximate 
length of the curve y = e* between x = 0 and x = 1. 


[Solutions to Exercises 1~6 on pp. 53-5] 
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Appendix: Solutions to the exercises 


Solutions to the exercises in Section 1 

1. @) Differentiating term by term, we have 
P'(x) = 8x? + 9x? + 2x — 7. 

(ii) Integrating term by term gives 


ee Sor Wa! 
Jocsrts = =x +9 t5* ge t2x+e 


where C is the constant of integration. 
Gi) Using nested multiplication, we have 
Up =2 
Uy = 0.29 + 3 = 3.58 
uy = 0.29u, + 1 = 2.0382 
uy = 0.29u, — 7 = — 6.408922 
Ug = 0.29uy + 2 = 0.141 41262. 
Hence p(0.29) = 0.141 41262. 
(iv) From (i), p(x) = 8x? + 9x? + 2x — 7, 
Using nested multiplication for this cubic polynomial gives 
Uy =8 
uy = 0.77up + 9 = 15.16 
Uz = 0.77u, + 2 = 13.6732 
us = 0.77uz — 7 = 3.528 364. 
Hence p'(0.77) = 3.528 364. 
2. Since f(x) = (1 + 8, we have 
Sey = Hl + xy? 
$x) = G)(-HL + x9? = 40 
When x = 0, (1 + x}* = 1 for any k. Thus 
£0) =1,f'(0) = and f"(0) = -4. 
Hence the 2nd Taylor polynomial is 
pix) = 1+ 4x — be? 
At x = 0.44 we evaluate p(x) using nested multiplication as 
Mo = —0,125 3 
uy = 0.44uo + 0.5 = 0.445 
uy = 0.44u, + 1 = 1.1958, 
ie. p(0.44) = 1.1958, 
This is correct to two decimal places, since 
(1 + 0.44)! = 1.2, 
3. f(x) =log,(1 + x). Differentiating gives 


y u ieee, 
S@=Ty "x)= aan" 
meee wy) a (yp! 
I"@)=tyap fe) = (Ie 
At x = 0 we have 
fO)=log.1=0, f0)=1, f"0)=—-1, 
FO) =2,-.-5 £0) = (— In — Yt 
Hence the nth Taylor polynomial about 0 for log.x is 
(-1)x? 2x3 (-197! (a — 1) 
PGI ees (oer apt nl 
+ 
=x-b? +$3-- eyes 


n 
as required. 


4. (i) Since 0 < c < 0.5, we have 1 <1 +¢< 1.5, and hence 
[l+el< 15. 

i) f(x) = (1 + x). Thus 
Se) = SL +x), "(x)= 201 + x), 
Sx) = 60(1 + xP? and f(x) = 120(1 + x). 

The error function for the 3rd Taylor polynomial is given by 


ox) = —¥ 646) —Sx4(1 + cx) 


where c, lies between 0 and x. 


Now we are given that 0 < x < 0.5, so it follows that 
0 <c, <0. and hence (by part (i)) that |1 + cx] < 1.5. Using 
this result to calculate a bound for e(x), we obtain 


Je(x)| < Sx* x 15=7.5x* (0<x<0.5) 
which is the required result, 
5. (i) We know that —1 <sin c <1 for all c, and hence 
{sin c| < 1. 
(i) Let f(x) = sin x. Differentiating gives 
Se) =cosx, "(x)= —sinx, 
fx) = —cosx, f(x) =sin x, 
SOx) =cosx and f(x) = —sinx. 
At x = 0, sin x = 0 and cosx = 1, Hence 
FO =0, SM=1, f"0)=0, 
f"O)=-1, £0) =0 and f°) =1. 
Hence the Sth Taylor polynomial about 0 for sin x is 
roelegece 
Px)=x- 3+ 5p 


The error is given using Taylor's Theorem, as 


6 
a(x) = Fines where c, lies between 0 and x. 


From part (i) we know that |sin c,| < 1. It follows that 
xs 
le) < 
Thus, at x = 1, we have 
le) < 4 = 0,001 388 89. 


Our error bound for the Sth Taylor polynomial about 0 
approximating sinx at x = 1 is thus 0.001 39, to three 
significant figures. (There is no point in quoting errors and 
error bounds to a large number of significant figures.) 

The actual value of sin 1 is 0.841 47099, and of the Taylor 
polynomial is 0.841 666 67, so the actual error is only 
0,000 196 (again to three figures). 


Solutions to the exercises in Section 2 


1. From the graph the three roots are approximately —0.5, 
—3.5 and —9.5. Using these values as initial estimates for 
the Newton-Raphson method with the recurrence relation 

_ (+1356 + 40)x, + 16.67 


or Gx, + 27x, + 40 
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we have the following table of results. (Note that the 
recurrence relation has been written in a form suitable for 
nested multiplication.) 


X, (Xo = —0.5) 


0} -05 =95 

1 | —0.49706422 —9.4570507 
2 | —0.497068 01 —9.456 5298 
3 | —0.49706801 —9.456 5297 


As you can see, after only 3 iterations the results agree to 7 
or 8 significant figures. Hence, to 5 significant figures, the 
roots are —0.49707, —3.5464 and —9,4565. 


2, (i) The equation 
xe* — 0,1 =0 


can be multiplied by 10e*, without introducing any new 
roots, to give 


10x — e* = 0. 


a 


Y¥' =3¥, + sinx, 


St 


3. The recurrence relation for the Taylor series method of 
order 1 is given by 


Wer =%+ hh. 
This is Euler's method (Unit 2, Subsection 2.2). 


4. The recurrence relation for the Taylor series method of 
order 2 is given by 


Wer =H +hh +4 ye 
Now the differential equation is 

y =3y+sinx. 
Differentiating gives 

y” = 3y' + cosx. 
Thus we have 

¥, = 3Y, + sinx,, 

¥," = 3¥/ + cosx,. 


¥," = 3%) + cosx, Wer = thy + hry? 


0 
0.02 
0.106 855 36 


0 
0.258 669 33 


Thus the original equation has roots whenever 10x = e*, 
Hence if we sketch the graphs of y = 10x and y = e* the 
intersections of the two graphs will give the roots (see Figure 
1). 


y 


Figure 1. Sketch of y = 10x and y = e*. 


There are two roots at approximately 0.2 and 3.5. 


(ii) The Newton-Raphson method for this problem uses the 
recurrence relation 


x, exp (—x,) — 0.1 


wt =~ exp (—x,))( =) 
ie, 
eee x,-0.1 SxP (x) 
1-x, 
r X,(Xo X,(Xq = 3.5) 
0} 02 35 
1 0.102 675 34 3.575 3819 
2 | 011174432 3.577 1511 
3 | 0.11183255 3.577 1521 
4 | 011183256 3.577 1521 


To 6 decimal places the roots are 0.111 833 and 3,577 152. 


0.02 


1 
1.756.0746 0.106 855 36 


The approximate solutions at x, = 0.2 and x, = 0.4 are 0.02 
and 0.106 855 36 respectively. 


5. The recurrence relation for the Taylor series method of 
order 3 is given by 


Yer =¥+h¥ +4 yy’ +hhy”, 
The differential equation is 


y =siny. 
Differentiating twice gives 
y" = y'cosy, 
y” = y" cosy — y’sin y. 
Hence the required equations are 
Yi, =Y+h¥ +40 ¥ +h? ¥,", where 
H=1 
¥ =sin ¥, 
y" = ¥/cos ¥, 


X" = ¥," cos ¥, — ¥,?sin ¥,. 
For x» = 0, Yo = 1 and h = 0.1 we have 
Yo = sin Yo = 0.841 47099 
Yo’ = Yocos Yo = 0.841 47099 x cos 1 = 0.454.648 72 
Y@" = Yo" cos Yo — Yo? sin Yo 
= (0.454 648 72 x cos 1) — (0.841 47099? x sin 1) 
= —0.350175 49. 


Thus 


Yi =1 + (0.1 x 0.841 47099) + (0.005 x 0.454 648 72) 
— (& x 0.001 x 0.350175 49) 


= 1.086 362. 


Solutions to the exercises in Section 3 


ope 


1. @ coshx= ~ Thus, 


ee 
© = -sinhx. 


sinh (—x) = 


‘The corresponding results for sines and cosines are identical, 
ie. 


cos(—x)=cosx and  sin(—x) = —sinx. 


d a eer) e&-eF 

(ii) qo az 2 } z = sinh x, 
Lee d [@&—-e"*\|_&+e* 
Aim = 2 z )- = cosh x. 


2 
The corresponding results for sines and cosines are almost 
identical, ie. 

Gia) moe x and f (eax) = —sinx 

dx dx é 
2. Taylor's Theorem for f(x) = sinh x gives the 2nth Taylor 
expansion about 0, with remainder, as 

3 
sinh x = sinh0 + xcosh0 + * sinh + cosh + -- 
3 xaatt 
+ amie? * Gre 

where c, lies between 0 and x. (Use the result of Exercise 1 


to do the differentiation.) With sinh 0 = 0 and cosh 0 = 1 we 
have 


‘cosh c, 


“ x me vals 
sinhxext+ stat + Gao 
sant 
aro 


As n becomes large we obtain the required expansion for 
sinh x, assuming (as is in fact true) that for fixed x the 
remainder term approaches zero for large n. 


3. y = Leosh ax. Thus, 


Ape 
anh 


‘Using the formula for the length of a curve, given in the box 
on p, 31: 


‘ re 
=f fis (2) ax 
: 
-{ J+ sinh? ax dx 
as 


a 
= |» cosh axdx 


gaps 


=} (in a —sinh (—a)) 


(since cosh? ax — sinh? ax = 1) 


=? sinha (using the result of Exercise 1). 
4, Substituting L = 2.08 into Equation (7) gives 


208 = 2sinh a 
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which can be rearranged as 
sinh a — 1.04a = 0. 


y=suha -l-Oha 


Figure 2. Sketch of y = sinha — 1.04a. 


From the sketch graph of y = sinha — 1.04a in Figure 2, we 
deduce that there is a root at approximately 0.49, The Newton- 
Raphson method gives 

pe Se sinh a, — 1.04, 

rei = a ~ V cosh a, — 1.04 | 
dy = 0.49 

a, = 0.48703033 

a; = 0.48700263 = as, 


Hence we have a = 0.48700263. 
From Equation (6) we have 
= 1 (cosha =-1) 
a 
= 0.248 35216, 
Thus the displacement for L = 2.08 m is, to 6 significant figures: 
d =~ 0.248 352m. 
5. The Sth Taylor polynomial for sinh a about 0 is 
a a 
P(x)=at etn 


Thus 
re 2 fai ce 
Lacanhass a+ 2+) 
=e 
3° 60 


Now to obtain a in terms of L in this approximation, we note 
that 
ea 
L=2+5+% 


is a quadratic in a’, ie. 


The formula method gives 
1 pt 
ged V9" @ 
1/30 
= 104/004 @L — 120 
= 10+ /ooL—20. 
As a? must be positive there is only one possibility, ic, 
@ = -10+2/15SL—5 
giving 


a=,/-10+2/15L—5. 


(L-2) 
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6, IfL = 2.08m we have from Exercise 5, remembering that it 
is an approximation: 


a ~ 048701859. 
Equation (6) gives 
d=+(cosh a—1) 
a 
~ 0.248 36061. 


Comparing this with the accurate result in Exercise 4, we see that 
the approximation is correct to 4 significant figures, and even the 
fifth significant figure is nearly correct. 


Solutions to the exercises in Section 4 
1. Let the equation of the quadratic be 
P(X) = do + ayx + ayx*, 
Since p(x) must pass through the three data points, we must have 
P(—1) =a — a, + a, =2 
PO) =a =3 
PU) =a9 +a, +a;=0. 


We can see that the solution is ag = 3,a, = —l anda, = —2,s0 
that 


P(x) = 3 — x — 2x7. 
2. The interpolating polynomial is of the form 
P(x) = ay + a,x. 


The closest data points to x = 1.134 are x = 1.13 and x = 1.14. 
Thus we require p(x) to satisfy the equations 


P(L13) = dg + 1.134, = 0.1222 

P(1.14) = dg + 1.14a, = 0.1310. 
This gives a, = 0.88 and a) = —0.8722, Hence 

p(x) = —0.8722 + 0.88x 
and 

p(1.134) = 0.12572. 
3. @) The maximum value of |e| in the interval —1 <c <1 
occurs at c = | since e* is a positive, increasing function in this 
interval, Hence in this interval, 

0<|e|<e=2.7183 to 4 decimal places. 
(ii) The equation for the error is given by 


e(x) = —(X — Xo)(X — x1)(x — ee 


where c, lies between xo and x3. 

Now f(x) = e*. Therefore f(x) = e*. Hence, 
efx) = —(x + Ix — SPs, 

Now for —1 <c, < 1 we have |exp(c,)| < 2.7183, and so 
fe(x)] < |(x + 1)x(x — ye 

For x = 0.5 this gives 


{e(0.5)| < [1.5 x 0.5 x (-05)748 


=0.170 to 3 decimal places. 


Thus an error bound for the interpolating polynomial 
approximation at x = 0.5 is 0.170. 


4. Euler's method in composite form for the integral of e* 
between x = 0 and x = 1 with h= 0.2 gives 


f "tx ~02(C9 +6? +04 4608 4 <8) 
F = 1.5522 to 4 decimal places. 
(The true solution is e — 1 = 1.7183 to 4 decimal places.) 
5. The composite trapezoidal method for the integral of e* 
between x = 0 and x = 1 with h = 0.2 gives 
Ih edx ~ 0.2 (he? +e? +, 04 + 0% + 0% + fel) 
3 = 1.7240 _ to 4 decimal places. 


The true solution is e — 1 = 1.7183 to 4 decimal places. The 
result is much more accurate than the result obtained using the 
composite Euler’s method. 


6. (i) To use Simpson's method to approximate the integral of 
e* between x = 0 and x =1 we have h = 0.5, giving 
7 
fede = Deer 4 409 401) 
0 
= 1.7189 to 4 decimal places, 


This is much more accurate, for much less work, than the 
approximations in Exercises 4 and 5. 


(i) With h=4 we have 
: 0.25 0.23 aE) 1 
ff ete (et se +205 + 40-78 4 el) 
o 
= 17183188, 


(The true solution is. 1.718 281 8.) 
With h = } we have 
[oe = M5 easing 200m 44ers 4 20s 
+ e525 4. 260-75 4 40-875 4 ol) 
= 1.7182842 to 7 decimal places. 
This is an extremely accurate result, 
7. The approximate value for 


2 
L= f (1 + cos? x)!/? dx 
lo 
using the composite Simpson's method with h = i is 


Lx = Ai + costo)! +4(1 + cos?) ” 
ane 3 


=" 


+2(1 + cos?= er 1 + cos’ 
4 8 


42 
2m 
+ ( + cos ;) } 
= 1.910141 to 6 decimal places. 


(A more accurate answer, obtained by computer, using 
Simpson's method with h = 7/200, is 1.910099 to 6 decimal 
places. The solution obtained is thus correct to 4 decimal 
places.) 


Solutions to the exercises in Section 5 
1. @) From the hint we know that 


2 tan x sec? x = 2 tan x (1 + tan? x) 
= 2 tan x + 2tan? x. 


34 


Since tan x is positive and increasing in the interval 
0 <x <0.1, the maximum value of the above function occurs 
at x = 0.1. Hence, 

2 tan c sec? c <2 tan 0.1 sec? 0.1 

=0.20 to 2 decimal places. 

(i) The linear Taylor polynomial about 0 is given by 

P(x) =f(0) + f'(0). 
For f(x) = tan x we have 


S') = sex = 


so that f(0) = 0 and f'(0) = 1. 


Hence the linear Taylor polynomial about 0 approximating 
tan x is given by 


P(x) = x. 
The error function is given by 
a(x) = —4x7f"(cx) 
where c, lies between 0 and x. 


Now f"(s) = sec! x= 2 tan x ee? x 


Hence 0 <f"(c,) = 2 tan ¢, sec? ¢y 
< 0.20 by part (i), since O< c, < 0.1. 
A bound for e(x) is therefore 
le(x)| < $x? x 0.20 = 0.1x?. 


Y¥=x,¥%+1 
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Hence for this problem the recurrence relation is 


Xero = Xe — (1) 


exp(x,)—3 


To obtain an initial approximation we sketch the graphs of 
y =e and y = 3x as in Figure 3. 


From the graph we can see that two curves intersect at 
approximately x = 0.6 and x = 1.5. With xo = 0.6, Equation 
(1) above gives 

xX, = 0.618778 47 

xz = 0.619061 22 

X3 = 0.619061 29 


=X 


With xo = 1.5, Equation (1) gives 


x, = 15123581 
xz = 15121346 
=X5 


Hence to 5 decimal places the two roots are 0.61906 and 
1.51213, 


3. The Taylor series method of order 2 uses the recurrence 
relation 


Yer =¥% +hY + 4h? ¥,", where 
¥=x¥+1 and 
Y=¥+x¥. 
With h = 0.1 we obtain the following table of results. 


¥+1=¥, + 0.1%" + 0.005Y," 


2. The Newton-Raphson method uses the recurrence 
relation 


Xe =~, fe. 


(x,) 
We have 
S(x)=e-3x, fix)=e-3. 
cn 
s+ 


— = 
O 02 OF 0608 HW RK MH x 
Figure 3. Graphs of y = e* and y = 3x. 


1.105 
1.222 1303 


Thus: when x = 0.1, y = 1.105, 
when x = 0.2, y = 1.222. 


4. (i) Let 
P(x) =a + bx 
be the interpolating polynomial. 


Since the graph of y = p(x) passes through the points (1,1) 
and (24, 14), we have 


pl) =a+b=1 

p(2k) = a + 24 = 1}. 
Solving for a and b gives a =} and b = 3. Hence 

p(x) = 06 + 04x. 
At x = 1.44 we have 

p(1.44) = 1.176. 
(ii) Since }x~*? is a positive decreasing function on the 
interval 1 < x < 24, its maximum value is at x = 1, and is 
equal to 4. Hence 

He3)<}  forl<c<2b 
(ii) From Theorem 4.1 on page 39, the error function is 
given by 

a(x) = (x — Ix — HE, 
Now f(x) = ,/x, and so 


S(xy=px"? and f(x) = —4x797. 
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From part (ii) we know that an upper bound for |f"(c.)| 
where I <c, < 2} is given by M = 4. Hence at x = 1.44, 


(1.44 — 1)(1.44 — 2.25), 
8 


fe(1.44)| < = 0.04455. 


The actual error at x = 1.44 is 
(1.44) = p(1.44) — \/1.44 = 1.176 — 1.2 
= —0.024 


and its magnitude is less than the error bound for the 
approximation. 


5. Simpson's method in composite form with h = 0.25, 
approximating the integral of tan x between 0 and 1, gives 
1 
f tanxdx = 925 (tan 0 + 4 tan 025 +2 tan 0 
0 


+ 4tan0.75 + tan 1) 
= 0,616 48052. 


(The true solution is 
1 
it tan x dx = [log. (sec x)]} = 0.615 62647.) 
lo 


6. From Section 3 the length L of a curve y is given by 


dy 
For y = 2. 
‘ory ¢* we have e& and 
1 
L=[ 1 + e** dx. 
0 
Simpson's method for this integral, with h = 0.25, gives 
pa 9B isa sire + T+e+4/i+e? 
+ ./1 + e?) = 20035275. 
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